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Abstract 

Developments in dynamical systems theory provides new support 
for the macroscale modelling of pdes and other microscale systems 
^T^ such as Lattice Boltzmann, Monte Carlo or Molecular Dynamics simu- 
• ^ lators. By systematically resolving subgrid microscale dynamics the 
^ dynamical systems approach constructs accurate closures of macroscale 
^ discretisations of the microscale system. Here we specifically explore 
reaction-diffusion problems in two spatial dimensions as a prototype of 
generic systems in multiple dimensions. Our approach unifies into one 
the modelling of systems by a type of finite elements, and the 'equation 
free' macroscale modelling of microscale simulators efficiently executing 
only on small patches of the spatial domain. Centre manifold theory 
ensures that a closed model exist on the macroscale grid, is emergent, 
and is systematically approximated. Dividing space either into over- 
lapping finite elements or into spatially separated small patches, the 
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specially crafted inter-element/patch coupling also ensures that the 
constructed discretisations are consistent with the microscale system/ 
PDE to as high an order as desired. Computer algebra handles the 
considerable algebraic details as seen in the specific application to the 
Ginzburg-Landau pde. However, higher order models in multiple di- 
mensions require a mixed numerical and algebraic approach that is also 
developed. The modelling here may be straightforwardly adapted to a 
wide class of reaction-diffusion pdes and lattice equations in multiple 
space dimensions. When applied to patches of microscopic simulations 
our coupling conditions promise efficient macroscale simulation. 

Keywords multiscale computation, discrete modelling closure, multiple 
dimensions, gap tooth method, centre manifolds, reaction-diffusion equations, 
computer algebra. 
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1 Introduction 

Computational simulation is a key enabling technology in engineering, science 
and other quantitative fields P, e.g.]. Coherent spatio-temporal dynamics 
is the preeminent example of complex system behaviour as it emerges from 
the interactions of many similar components at each locale in space [151 EH 
[2T| \T6\ e.g.]. We must simulate such systems on the scale of interest and 
operation. But systems that depend on physical processes over multiple 
scales pose notorious difficulties. These multiscale difficulties are major 
obstacles to progress in fields as diverse as environmental and geosciences, 
climate, materials, combustion, high energy density physics, fusion, bioscience, 
chemistry, power grids and information networks [8]. Following the 'equation 
free' approach of Kevrekidis, Samaey and colleagues [2Sl e.g.], we here address 
the extraction, using dynamical systems theory, of computationally efficient 
macroscale models from given microscopic models, whether pde or lattice 
dynamics or other microscale simulators. 

The ultimate aim of this article is to provide theoretical support for and to 
further develop Kevrekidis' et al. 'equation free' approach to multiscale 

modelling. Given a numerical simulator for physical components at much 
smaller scales than the scale of primary interest, the 'equation free' aim is to 
bridge space and time scales to simulations resolving the macroscale of primary 
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interest. Here we bridge space scales by generalising to multiple dimensions 
(specifically 2D) both the methodology and supporting theory for the 'equation 
free', gap-tooth method for microsimulators that was initiated by Gear, Li 
& Kevrekidis and Samaey, Kevrekidis & Roose [51], |52] ■ Figures [T]-[3] 
show snapshots of an example simulation of the gap-tooth method in 2D: 
microscale simulators executing on a fine grid only within the 64 patches 
(of an 8 X 8 macroscale grid) are coupled across empty space to efficiently 
simulate the dynamics over large spatial scales; the computational cost here 
is one sixth that of a microscale simulation over the whole domain, but 
much better gains may be obtained (but do not provide suitable graphics). 
Crucially, although our analysis considers systems in principle representable 
in the class of the reaction-diffusion pde ([T]), in application the gap-tooth 
method does not require knowledge of the specific PDE: the gap-tooth method 
provides on the fly closure. Such closure constitutes critical components of, for 
example, mathematical homogenization |52 ] [T7 ] IT| e.g.], renormalization group 
techniques (T0l|37lini e.g.], and multiscale finite elements [201 H] e-g-]- But by 
avoiding the need to algebraically find the closure the gap-tooth scheme has 
the potential to efficiently simulate many systems whose macroscale dynamics 
are otherwise unaccessible. 

The key to support the gap-tooth scheme of Kevrekidis et al. ||24j is to couple 
patches of spatial dynamics across space; Figures [l]-[3] show an example. Such 
coupling needs to preserve the accuracy and stability characteristics of the 
microscale dynamics. Here we use the dynamical systems theory of centre 
manifolds [S] E] 155] e.g.] to guarantee a controlable level of fideUty between 
the microscale and the macroscale simulation. Thus the first contribution 
of this article is to extend the dynamical systems approach of the so-called 
'holistic discretisation' e.g.] from one spatial dimension to the discrete 

modelling of the class of two dimensional, homogeneous, nonlinear reaction- 
diffusion equations 

-^ = V-[f(u,Vu)Vu]+ag(u), (1) 
for the field u(x,ij,t). In principle we could consider the pde on any typical 
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u(x,y,t) at t=0.0000 




x/H 




Figure 1: perspective view of 8 x 8 
patches of a microsimulator on a 
macroscale grid of spacing H = 9.6 
at nondimensional time t = . 
Each patch is a microscale discreti- 
sation of the Ginzburg-Landau 
PDE ([2]) with nonhnearity a = 3 
on a 5 X 5 fine grid: the microscale 
simulator executes on only 1 6% of 
space. This initial condition has 
random significant microscale fluc- 
tuations superimposed upon a large 
scale variation. This initial condi- 
tion evolves to Figures [2] and [3j 
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u(x,y,t) at t=0.2500 





x/H 



Figure 2: perspective view of 8 x 8 
patches of a microsimulator on a 
macroscale grid of spacing H = 9.6 
at nondimensional time t = 0.25 . 
Each patch is a microscale discreti- 
sation of the Ginzburg-Landau 
PDE ^ with nonhnearity a = 3 
on a 5 X 5 fine grid: the microscale 
simulator executes on only 16% 
of space. At this time the spa- 
tial fluctuations within each patch 
have nearly smoothed to reflect the 
macroscale variations. 
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u(x,y,t) at t=1.0000 




x/H 



(view u3d in Adobe Reader) 



Figure 3: perspective view of 8 x 8 
patches of a microsimulator on a 
macroscale grid of spacing H = 9.6 
at nondimensional time t = 1 . 
Each patch is a microscale discreti- 
sation of the Ginzburg-Landau 
PDE ([2]) with nonhnearity a = 3 
on a 5 X 5 fine grid: the microscale 
simulator executes on only 16% 
of space. By now the variations 
within each patch are smooth and 
the patch evolution reflects the dy- 
namics of macroscale pattern for- 
mation. 
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domain with Dirichlet, Neumann or mixed boundary conditions |l3]; however, 
for simphcity, in this article we generally restrict attention to spatially periodic 
solutions so that the modelling is homogeneous in space. Generalisation to 
spatial dimensions higher than two appears straightforward. 

To achieve this aim of providing effective systematic closures for macroscale 
models in multiple dimensions, our approach systematically models subgrid 
microscale processes. For example, continuing the gap-tooth simulation of 
Figures [l]-[3] would enable reasonable exploration of the competition between 
meta-stable macroscale domains where u ±1 in the Ginzberg-Landau 
PDE ([2]). Sections [2] and [6] discuss two distinct avenues of theoretical support 
for our modelling: respectively that of centre manifold theory [HI EH I2SI El 
e.g.], and that of classic consistency. Such dual justification is a strength of 
this approach. 

The complex dynamics we address arise through the interaction over space of 
local microscale dynamics whether in a PDE such as ([T]), or a discrete lattice 
equation [211 [151 e.g.], or a microscale simulator (Figures [I]-[3]). Analysing the 
dynamics of a PDE for fixed macroscopic grid spacing, a third contribution of 
this article is to use centre manifold techniques to underpin accurate models 
of nonlinear dynamics by resolving naturally the dominant subgrid microscale 
structures and their interactions, both internally and with macroscale vari- 
ations. Instead of imposing a subgrid field, such as the usual polynomial 
interpolation of finite differences and finite elements, here the PDE deter- 
mines the subgrid field. Then the derived macroscale closures enable a 
relatively coarse numerical grid to significantly improve computational speed 
and stability in numerical solutions of the PDE. For example, we expect more 
extreme parameter regimes may be explored without the need for artificial 
hyper- viscosities [53 e.g.]. 

The analysis of pdes in Sections [2]-^ parallels, and has much commonality 
with, our analysis of dynamics on disjoint spatial patches in Sections |2] 
and [5j A micro-simulator within each patch requires boundary values. If 
the microsimulator were to be executed over the entire macro-domain, then 
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such boundary values come naturally from immediately neighbouring fine 
grid points; such neighbours are missing in gap-tooth simulation such as 
those of Figures [l]-[3} Instead we propose the innovation that classic Lagrange 
interpolation from surrounding macroscale grid values provides the accurate 
coupling for the small microscale patches, analogous to accurate coupling of 
one dimensional dynamics [HI HQ]. The centre manifold theory of Section [2] 
then supports the macroscale modelling. To complement this dynamical 
systems support, Theorem [6] provides support for the consistency of the 
approach: the order of consistency growing linearly with the order of the 
interpolation. This classic coupling rule establishes a strong connection 
between classic finite difference discretisations of pdes, finite elements, and 
the methodology of the gap-tooth method. 

Our work here presents two faces to computational simulation. On the one 
hand we present a preprocessing methodology for generating potentially highly 
accurate closures of discretisations of pdes or lattice dynamics. These would 
closures would subsequently be used to markedly speed up simulations. On 
the other hand we prove that the same approach provides coupling conditions 
for accurate and effective, on the fly, closures for the 'equation-free' macroscale 
simulation of highly detailed microscale dynamics]^ 

As a particular example, previewed in Figures [T]-[3[ Sections [3[ |4] and [5] explore 
in some detail the modelling of the real valued, two dimensional, Ginzburg- 
Landau equation obtained from the pde ([T]) with cubic reaction, g = u — u^, 
and constant diffusion, f = 1 , namely 

^ = V^u+ a(u-u^) . (2) 
9t 

We choose this 2D real Ginzburg-Landau equation as an example prototype 
reaction-diffusion pde because its dynamics are well understood [HI 123 e.g.]. 



^In either case, the issue of parallehsing the computational simulations are the same, 
and familiar from usual approximations: to obtain higher order accuracy, generally use a 
wider computational stencil, which requires proportionally more communication between 
parallel processors in some domain decomposition of the computation. 
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and because this pde is important as a phenomenological model [291 P-6, 
e.g.]: much interest hes in the long time evolution of the interacting domains 



of the quasi-stable states u ~ ±1 . In Section [4.2| the steady states of the 
example2D Ginzburg-Landau equation ([2]), and their stability, measure the 
accuracy and effectiveness of various order models in this application. Sec- 



tion 3.1 briefly compares a low order model with a classic flnite difference 
model to indicate their similar performance and to provide a base for compar- 
ing high order models. MacKenzie [33] reported that our modelling, based 
upon a mosaic of local dynamics, is much less costly to use than a global 
description of an inertial manifold [531 ESI |22l e.g.]. 

The macroscale model is based upon dividing the domain into either over- 
lapping flnite elements or disjoint small patches. Following analogues in 
one dimension [SI |32l e.g.], neighbouring elements/patches are coupled with 
strength parametrised by y. Section [2721 then discusses how centre manifold 
theory [3l [551 [261 [3 e.g.] assures us of the existence of an exactly closed dis- 
crete model. Further, this discretisation is exponentially quickly attractive — it 
contains the emergent dynamics. Although we cannot exactly construct this 
closure, theory asserts it may be approximated to any order in the strength 
of the inter-element/patch coupling y and the nonlinearity a. 

The special coupling conditions Q assure us that the resultant macroscale 
discrete models are also consistent with the dynamics of the reaction-diffusion 
PDE (Section [g]). A further contribution of this article is that the proofs for 
consistency in Section [6| are new and more powerful, leading to new theorems 
on nonlinear and multidimensional consistency, and to a new theorem on the 
consistency of patch dynamics. 

Section [3] outlines the construction, consistency and predictive accuracy of 
low order asymptotic approximations to the macroscale discretisation of the 
Ginzburg-Landau pde ([2]). To extract another order of accuracy from the 
algebra, we flnd (for the flrst time) the adjoint operator of the diffusion 
operator on the elements/patches with the nonlocal coupling conditions. The 
null space of this adjoint, strikingly similar to a Galerkin basis, enables us 
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to use an integral solvability condition to construct the third order discrete 
model. 

However, we cannot algebraically find higher order models nor any of the 
patch models. This inability to construct algebraic approximations is one 
major difference between systems in one and multiple spatial dimensions. In 
a further contribution, Section |4] introduces how to numerically construct 
the microscale subgrid field and its evolution in 2D reaction-diffusion pdes 
using the Ginzburg-Landau pde as an example. Such integration of numer- 
ical solutions for the subgrid field in complex algebraic expressions for the 
macroscale parametrisation is novel. We find that even a relatively coarse 
subgrid microscale resolution is adequate to reasonably accurately underpin 
the macroscale modelling. 

2 Divide the domain into elements/patches 

We place the discrete macroscale modelling of general, two dimensional, 
reaction-diffusion dynamics within the purview of centre manifold theory by 
dividing the domain into either overlapping square elements or into small 
disjoint separated patches, as shown schematically in Figure |4j 

2.1 Extend non-local coupling conditions to multiple 
dimensions 

Define a grid of points (X^, Yj) with, for simplicity, constant spacing H, see 
Figure |4] The i, jth element, Ey, is centred upon (Xi,Yj) and of width 
Ax = Ay — 2rH : when r = 1 we cater for the overlapping elements of holistic 
discretisation [44^ e.g.]; when r < 1 /2 we cater for the spatially separated 
patches of equation free modeUing jlHl 1251 e.g.]. Let xiij(x,'y,t) denote the 
field in the i, jth element /patch. The fields uy(x,i|,t) evolve according to 
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Figure 4: Discretise a 2D domain into square elements/patches. The i, jth ele- 
ment/patch (solid blue/magenta) is centred upon the grid point (Xi, Yj): when 
r = 1 (bhie) Eij overlaps neighbouring elements to extend to the neighbouring 
grid points; when r < 1 /2 (magenta) E^j forms a patch separated by empty 
space (gaps) from neighbouring patches. 
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the reaction-diffusion PDE ([T]); that is, 

9x4-' 

^ = V ■ [f (Ui^j, Vut^j) Vuy] + cxg(ug) . (3) 
The original field u(x,i|,t) is then predicted by itij(x,i|,t) for (x,!)) G Ey. 

The evolution of the field over the whole domain depends upon how the 
elements/patches are coupled together. To couple the dynamics of each 
overlapping element to its neighbours, the case r = 1 , we use 'inter-element 
coupling conditions' (iCCs) around the i, jth element of 

u-i,j(Xi±i,-y,t) =yiti±i,j(Xi±i,y,t) + (1 -y)uy (Xi,-y,t), [ij - Yj| < H, 
u-i,j(^,Yj±i,t) =yuy±i(x,Yj±i,t) + (1 -y)uy(x,Yj,t), |x-Xi| < H. 

(4) 

These iCCs are a natural extension to 2D of iCCs established for ID dynam- 
ics im e.g.]. The crucial feature is: with y — the elements are effectively 
isolated from each other, dividing the domain into decoupled elements with 
consequently independent dynamics; whereas with y = 1 these ICCs ensure 
sufficient continuity between elements to recover the original problem over 
all space. The ICCs (|4]) embeds the physical problem, parameter y = 1 , into 
a family of problems, general y, and then we access the physical problem 
from the tractable base at parameter y = . Modelling via these overlapping 
elements is called 'holistic' because within these elements we resolve subgrid 
structures by systematically approximating solutions of the PDE (|3]); that 
is, the PDE itself tells us what are appropriate subgrid fields. In contrast, 
methods such as finite differences, finite elements and finite volumes, im- 
pose an assumed subgrid field upon the elements (typically a relatively low 
order multinomial). Interestingly, the coupling Q of overlapping elements 
appears to have analogues in other multiscale methods: the 'border regions' 
of the heterogeneous multiscale method (Qj e.g.], the 'buffers' of the gap-tooth 
method [521 6.g.], and the overlapping domain decomposition that improves 
convergence in waveform relaxation of parabolic pdes [121 e.g.]. 

We now consider the multiscale case of how to couple patches across space, 
as shown in Figures [l]-[3j To couple the dynamics of each separated patch 
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to its neighbours, across empty space as this is the case r < 1 /2 , we use the 
different couphng conditions (iCCs) around the 1, jth patch of 



uy (X, ± rH, y , t) = Sf' (y )£:f ^ (y )u,,j (X^ , Yj , t) , j-y - Yj | < rH , 



uy (x, Yj ± rH, t) = (y )£:f '■(y )uy (X^, Yj , t) , |x - XJ < rH 



(5) 



in terms of subgrid variables £, = (x — X^j/H and "H = (y — Yj)/H, and the 
discrete shift operator £. The reason for the difference between the iCCs (|4]) 
and (|5| arises because the patch iccs ^ depend only upon the neighbouring 
grid values Uy(t) = Ut j (Xi, Yj, t): in large scale computational simulation, 
the ICCs ([5]) minimise communication between patches in comparison with the 
ICCs Q which require data along the mid-patch lines to be communicated with 
neighbouring elements. The shift operators <5^^(y) and E^'^ly] on the right- 
hand side of the ICCs ([s]) derive from classic Lagrange interpolation through 
neighbouring grid values, but ameliorated by the coupling parameter y. These 
y ameliorated shift operators are best defined in terms of the well known 
centred difference 5 and mean \i. operators: for each direction indicated in ^ 
by 1 and j, the y ameliorated shift by an fraction ±£, of a grid spacing is 
represented in classic operator algebra [38l [19], e.g.] as 

= 1 + £,y(±^5 + \b^) + 1 )y'(±^5 + \b^Y + ■ ■ ■ . (6) 

Identically for the ICCs (|4]), a crucial feature of the ICCs ([s]) is: with y — 
the patches are effectively isolated from each other, dividing the domain into 
decoupled patches with consequently independent dynamics; whereas with 
y = 1 these ICCs ensure each patch matches those in its neighbourhood via 
classic Lagrange interpolation. 

Whether patches or overlapping elements we need domain boundary conditions 
to form a well posed problem. For definiteness in theoretical support, let there 
be m elements/patches in the spatial domain £1 with the field required to be 
periodic in both x and y, and the field to be in a Sobolev space such as H^(0). 
For example, the elements/patches may form a -\/ttl x grid in the domain 
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(any factorisation of m is feasible). In principle, the 2D elements/patches 
could be any shapes, regular or irregular; square ones appear to be easiest 
to start with. Although physical domain boundary conditions have been 
explored for ID domains |l3l EU e.g.], in this initial work we avoid issues 
of physical domain boundary conditions on the domain Q by the doubly 
periodic condition. Such doubly periodic boundary conditions enable the 
most straightforward theoretical statements of support. 



2.2 Centre manifold theory supports discrete models 

This section describes in detail how the iCCs Q or (|5]) lead to centre manifold 
theory [31 ESI I2SI El e.g.] supporting an accurately closed, discrete model for 
general reaction-diffusion systems ([T]) via its dynamics (|3]) within coupled 
elements. Figure 11 shows one example of the scale separation that underlies 
this modelling: the figure shows a large spectral gap between the 'small' 
eigenvalues (decay rates) of the macroscale modes and the eigenvalues < —40 
of the rapidly decaying microscale modes. 

A homotopy in the coupling parameter y connects the physically relevant 
macroscale discretisation to a theoretically tractable base. When parameters 
a = y = both the pde reaction and the coupling on the right-hand side 
of the ICCs (|4]) or ([5]) disappear. The elements/patches are then effectively 
isolated from each other and so the resultant diffusion in the pde ([s]) is 
particularly simple: exponentially quickly in time, the field Utj becomes 
independently constant within each element. We use this m parameter family 
of piecewise constant solutions as a basis for analysing the case when the 
elements/patches are coupled together, y 7^ . Particularly interesting is 
the approximation for full coupling, y = 1 , when the pde ([T| is effectively 
restored: over the whole domain because ICCs (|4]) then ensure sufficient 
continuity between adjacent elements as described previously for pdes in one 
spatial dimension |l2l [321 HHl e.g.]; or because the ICCs (|5| connect patches 
over the gaps via clasic Lagrange interpolation. 
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The support of centre manifold theory j3l [551 ISSl O e.g.] is based upon a 
linear picture of the dynamics. Adjoin to the pde ([T| the dynamically trivial 
equations 

^ = ^ = 0, (7) 

9t at ' ^ ' 

and consider the reaction-diffusion dynamics in the extended state space 
(uijfxjijjjy, a). In this extended state space, points a = y = and 
TXij = constant, say Utj, are equilibria of the diffusion ([s]), hence these 
form a subspace of equilibria in the extended state space, the subspace 
Eo = {(Uij, 0, 0) I for all Uij}. Linearized about each of the equilibria in Eq, 
the PDE for perturbations u(j(x, ijjt) within each element /patch is then the 
constant coefficient diffusion pde 

= fyV u^j for (x,y) G Eg for each (8) 

where the constant diffusivities fij — f(Uij,0). These pdes are decoupled 
because they are to be solved with the y = iCCs 

u(j (X, ± rH, y, t) = u(^j (X^, {y, YJ, t), |xj - Yj| < rH , 
u(j(x,Yj±rH,t)=u(j({x,X^},Yj,t), |x-Xt|<rH, 

where the first alternative in braces on the two right-hand sides corresponds to 
the elemental iCCs (|4]), whereas the second alternative in braces corresponds to 
the patch iccs ([s]). Among other eigenmodes such as Ug = 1 with eigenvalue 
zero, separation of variables shows that the following are some of the linear 
eigenmodes associated with each element/patch: 

a = y = , u( j oc e'^^-'-'''^^ sin(k7r£,/r) sin(l7rri/r), 

inside the i, jth element for all integers k, I = 1 , 2, 3, . . ., where the eigenvalues 
corresponding to each of these modes are 

Ay ,k,i = -fg . (10) 
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Simple numerical solutions of the pde (jsj) and boundary conditions ([9 



confirm that ( 10 ) for k, I = 0, 1 , 2, 3, . . . are the only eigenvalues on elements 



(k, I = included here), but that on patches the eigenvalues, in addition to ( 10 ) 
for k,l = 1,2,3,..., are l%^7^/[r^H^) for % = 0,-1.250,-3.250,-3.73 ± 
12.00, These eigenvalues Afij7T^/(r^H^) of the decoupled dynamics ap- 



proximate the clusters shown by the histogram of Figure 11 for the fully 
coupled dynamics. In addition there are the two trivial 'parameter' eigen- 
modes: firstly y = constant and a = Ug = ; and secondly, cx = constant 
and y = u( j = . In a spatial domain with m elements / patches and when 
all diffusivities ftj > 0, then all eigenvalues are negative, — fij7r^/(r^H^) or 
less, except for m + 2 zero eigenvalues. Of the m + 2 zero eigenvalues, one is 
associated with each of the ra elements/patches and two come from the trivial 
equations ([T]) for the parameters. That is, the slow subspace is {(ud j,y, a,)} 
for constant Uij. The above spectrum establishes the following corollary of a 
centre manifold existence theorem [31 |55l p. 281, p. 96 respectively]. 

Corollary 1 (Existence) Provided the nonlinear diffusivity f and reaction g 
in ([3]) are sufficiently smooth, and all f^j > then an m + 2 dimensional 
slow manifold Ai exists for ^ and Q, or ([s]) and in some finite 
neighbourhood of the subspace Eq of equilibria^ 

The slow manifold A4 is parametrized both by the two parameters y and a, 
and by a measure of the field in each element /patch; we use the grid value 
Ui,j(t) = iXij(Xi, Yj, t) to measure the field in the i, jth element /patch. Using 
U to denote the vector of such parameters, we write the slow manifold Ai as 
the subgrid fields 

itij =Ui,j(x,-y;U,y, cx). (11) 

■^The eigenproblem for the pde and boundary conditions were solved via second order 
finite differences on microscale grids of size 9x9, 17x17 and 33 x 33. Then the eigenvalues 
were extrapolated with a Shanks transform to approximately the reported accuracy. 

^Keep clear the distinction between centre manifold theory and the slow manifolds 
discussed here: centre manifold theory applies to systems where the real part of the 
eigenvalues of critical modes are zero; whereas here we explore and construct the particular 
case of slow manifolds because here the relevant eigenvalues are precisely zero. 
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These functions xii jfx, 'y;U,y, a), that Sections [s] and [i] construct for the 
Ginzburg-Landau pde ([2]), resolve the subgrid scale physical structures as a 
function of the neighbouring grid values in U. On this slow manifold Ai the 
grid values Uij evolve deterministically: 

dUy/dt = U,,j = g,,j(U,y,a), (12) 

where gij is the restriction of ([s]) and Q, or (|3| and (|5|, to the slow mani- 
fold M.. In essence, this closure of the grid scale dynamics comes from the 



resolution of subgrid scale structures. The set of odes (12) form the discrete 
model of the original pde. 

Using the value of the field at the grid points to parametrise the slow manifold 
provides the necessary 'amplitude conditions' to close the problem: 

Ug=Ui,j(Xt,Yj;U,y,a). (13) 

Many other amplitude conditions are possible such as defining the 'ampli- 
tudes' Uij to be the mean field over the 1, jth element /patch. However, using 
the grid values are simple, traditional, and have a direct physical meaning. 

Centre manifold theorems |3l [551 6.g.] also support the following crucial emer- 
gence (called asymptotic completeness by Robinson [50j) and approximation 
properties. 

Corollary 2 (Emergence and approximation) Provided the nonlinear 
diffusivity f and reaction g in ^ are sufficiently smooth, and all fy > , 
then 

• every solution of the nonlinear reaction- diffusion dynamics ^ and 
or ^ and ([s]), that stays within a neighbourhood of the slow mani- 



fold Ai, (11), approaches exponentially quickly a solution of the discrete 



model (12) on the slow manifold (11); and 



the order of error in asymptotically approximating the slow manifold 



and its evolution, (11)-(12), is the same as the order of residuals of 
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the governing equations ^ and Q, or ([3]) and ([5|. In particular, 
because the base equilibria form a subspace Eq, here the approximation 
is global in the grid values Uij, it is local only in the coupling y and 
nonlinearity a. 



3 A slow manifold discretisation of the 
Ginzburg-Landau PDE 

We now explore constructing slow manifold, discrete models for a specific 
example reaction-diffusion system, the 2D Ginzburg-Landau PDE ([2]). In 
applications this construction is a preprocessing step prior to large scale 
numerical simulation. The support of centre manifold theory for our exotic 
closures suggests subsequent numerical simulation will be significantly more 
accurate and/or efficient. 

This section only addresses the slow manifold discretisation on overlapping 
elements with the aim of deriving and testing accurate nonlinear closures of 
the reaction-diffusion dynamics; later sections return to the slow manifold 
view of the dynamics on spatially separated patches. In contrast to one 
spatial dimension where it is straightforward to construct slow manifold 
discretisations for general pdes (see our web service US]), one obstacle 
in multiple dimensions is the limited range of slow manifolds representable 
by multi-variate polynomials. Consequently, the next section extends the 
approach via numerical solution of the subgrid fields. The subgrid fields on 
patches require such numerical solutions and so patches are left for later 
sections. 



Substituting the slow manifold ansatz ( 11 ) and ( 12 ) into the Ginzburg-Landau 



PDE (|2]), assuming solutions are doubly periodic, we obtain the PDE to solve 
for the model by equating the PDE for 9uy/9t to that obtained by the chain 
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rule: 

"^=11 |S^9Ki = V^u,, + a (u,, - uy . (14) 



k,l 



To construct the slow manifold ( 11 )-( 12 ) by solving the PDE ( 14 ) with coupling 



and amplitude conditions involves considerable algebraic detail. Computer 
algebra handles the algebraic details of the construction by iteration jUl e.g.]. 
The precise procedure used here is fully documented on a freely accessible, 
preprint server [M] to empower readers to check, reproduce and build upon 
our quoted results. The methodology constructs a model by driving to zero. 



to some order of error, the residuals of the governing differential equation (14) 
and the inter-element/patch coupling iccs Q or ([s]). By the Approximation 
Corollary |2] we construct correspondingly accurate approximations to the slow 



manifold of (14). These approximations, upon setting coupling parameter 
y = 1 and the nonlinearity parameter a appropriately, form 2D discrete 
models of the 2D Ginzburg-Landau PDE ([2]). 

Recall that although the approximations to the slow manifolds considered 
here are asymptotic, the existence of an emergent slow manifold holds in a 
finite domain around the slow subspace Eq (Corollaries [I]-[2]). Thus although 
we construct and label the slow manifold approximations via asymptotics 
in small parameters, the resulting models hold for finite parameter values. 
The issue is how well do the approximations perform: we give evidence that 
reasonable results are straightforwardly obtained for low order approximations 
for nonlinear parameter a up to about 30. 

One consequence of using computer algebra is that there is no need to record 
in this article most of the considerable algebraic detail in constructing the 
models. Those wishing to verify the correctness of the results recorded herein 
should download and examine our freely available technical report [M] that 
details the precise computer algebra procedure. Because the algorithm is 
based upon driving the residuals to zero, the critical aspect of the procedure is 
simply the correct coding of the computation of the residuals of the governing 
equations. One may straightforwardly edit the code j3l] to construct holistic 
discretisations for other react ion- diffusion systems in the class ([T|. 
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The 0{y^ + a^) holistic discretisation Satisfying the PDE and iCCs 
to residuals of 0{y^ + a^) the computer algebra procedure [Ml §2.2] gives 
subgrid fields which are too complex to record here. The corresponding 
evolution of the grid values on the slow manifold are 

+ [3(^ij5jUt,j)^ + l(6fUt,i)^] (2 + 6f) + (5fUt,j)^}Uy 
+ 0(y^ + a3), (15) 

where the operator difference 6 = iS^^^ — 8^"^^^ and mean \i = [E^l'^ + £^^^^]/2 
in each of the two i and j directions, and where the bold centred difference 
operator applies in both spatial dimensions: 

= U^+i,j + U,_i,j + Uy+i+Ug_i-4Uy, (16) 

= Ui+2,j + Ui_2,j + Uy+2 + Uy_2 

- 4(U^+i,j + U,_i,j + Uy+i + U,,j_i) + 12Uy . (17) 



The model (15) is simply the extension to two spatial dimensions of the 
0{y^ + a^) holistic model of the ID Ginzburg-Landau PDE [33]. 

The holistic discrete model has the dual justification of consistency with the 
PDE in addition to the justification provided by centre manifold theory. As 
proven in Section [6| consistency for such discrete models follows from the 
coupling ICCs ^ [44j. Set the coupling parameter y = 1 in the discrete 
equation (15) to recover the holistic discrete model of the Ginzburg-Landau 
PDE ([2]) in 2D. To test consistency, we expand the finite differences of (15) in a 
Taylor series in the grid spacing H [33', §3.2] to find the equivalent continuum 



PDE for the Oiy^ + (X^) holistic model (15) is 

VWa(u-u=) + ii^u|VuP-^ (1^ + g)+0(a^+H') . (18) 
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The 0{y^ + a^) holistic model is (9(H^ + cx^) consistent, maintaining in 2D 
the dual justification of holistic discretisation found for ID pdes HI]. 
Section ^ proves this consistency in some generality for 2D pdes (IT]). 



3.1 Illustration of the subgrid field in 2D 

Here we plot the subgrid fields for a coarse grid solution of the 0(y^, oc^) holis- 
tic model (obtained from (15) by omitting the term and then evaluating at 
y = 1) of the 2D Ginzburg-Landau PDE (|2|. We restrict attention to doubly 
odd symmetric solutions that are also doubly 27r-periodic: that is, 

( ^\ /-^f27r-x,ij,t) = -u(x,27r-y,t), odd symmetry, 
u(x,-L|,tJ = < (19) 
I u(x + 27t, y , t) = u(x, y + 27r, t) , 27t-periodicity. 



Figure [5] shows the subgrid fields for the 0(y^, cc^) holistic model with 4x4 
elements on [0, tt] x [0, 7t] at nonlinearity a = 6. The subgrid fields exhibit 
the nonlinear subgrid structure of the 2D Ginzburg-Landau PDE ([2]) and 
its interaction through the iCCs. The subgrid fields are comprised of actual 
solutions, albeit approximate, of the 2D Ginzburg-Landau PDE ([2]). 

Note the subgrid fields have small but noticeable jumps across the boundaries 
of the elements. Higher order holistic models reduce these jumps across the 
boundaries as seen, for example, in the holistic models of the ID Kuramoto- 
Sivashinsky equation 



The C>(y^, a^) holistic model needs improving We investigate the 
performance of the 0(y^, (xf) holistic model on coarse grids by comparing its 
bifurcation diagram to an accurate solution. Again we restrict our attention 



to doubly odd symmetric solutions that are 27r-doubly periodic (19). 



Continuation software auto [7| and xppaut [llj calculates this bifurcation 
information, as outlined for the Kuramoto-Sivashinsky equation in MacKen- 
zie's PhD dissertation |33]. In such bifurcation diagrams the blue curves 
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^ X 

Figure 5: an example of the subgrid field for the 0(y^, oc^) holistic model 
with 4x4 elements on [0, tt] x [0, tt] at nonlinearity a = 6 



Tony Roberts, January 13, 2013 



3 A slow manifold discretisation of the Ginzburg-Landau PDE 



24 



2D-RGL 0(/,a^) holistic, 8x8 elements 

8 1 ■ ■ 1 




a 



2D-RGL 0(7^, a^) holistic, 8x8 elements 

8 1 ■ ■ 1 




10 20 30 



a 



Figure 6: Bifurcation diagrams for the 2D Ginzburg-Landau PDE pi) with 
8x8 elements on [0, tt] x [0, 7t] for hoHstic models (a) C(y^, oc^) (b) C(y^, a^) . 
The accurate bifurcation diagram is shown in grey. 



indicate stable steady state solutions and red curves indicate unstable steady 
state solutions. The open squares indicate steady state bifurcations. 

Underlying Figure [6j in grey for reference, is an accurate bifurcation diagram 
of the 2D Ginzburg-Landau PDE (|2]). It is constructed with a computationally 
expensive, fourth order, centered difference, approximation with 24 x 24 points 
on [0,7r] X [0, tt]. The trivial solution undergoes steady state bifurcations 
at <x = 2,8,18 leading to the unimodal, bimodal and trimodal branches 
respectively. For nonlinearity ranging over 1 < a < 30, only the unimodal 
branch is stable; all other branches are unstable. 

Figure [6]( left) compares the bifurcation predictions of the 0(y^, oc^) holistic 
model with the exact bifurcation diagram (grey). The predictions are reason- 
able for parameter a up to about ten, and for the smaller range of amplitudes 
shown. Figure [oj^right) shows the significantly improved accuracy of the 
(9(y^,a^) holistic model that is obtained by resolving more inter-element 
interactions through retaining higher orders in the coupling parameter y. 
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As proven for other PDEs in ID [l2l |3T1 [32], higher order modeUing in 2D 
evidently improves predictions. 

Higher order models need numerical construction To improve the 
accuracy of such discrete closures we need to compute higher orders in either 
coupling y, or nonlinearity cx, preferably both. Improved accuracy occurs 
at higher order in comparable ID problems |32]. However, apparently it 
is not possible to analytically construct higher order subgrid fields in 2D: 
the subgrid fields required for our closures appear to be no longer in the 
class of multivariate polynomials. For example, the 0(y^, (X^) model used 
for Figure [6]^ right) is not immediately obtainable from the computer algebra 
iteration [31]. Instead, numerical methods must be used to find the subgrid 
fields as described in Section|4j Section [5] employs the same numerical methods 
to also construct slow manifold models on patches. However, the well known 
'solvability condition' in asymptotic mathematical methods empowers us to 
derive the next order in the evolution, analytically from the residuals, without 
needing to find the next order of the subgrid fields, and hence provides the 
model underlying Figure |6]( right). 

3.2 The adjoint provides an extra order of accuracy 

We scrounge an extra order of accuracy from the 'solvability condition' |39| 
e.g.] applied to residuals of the next asymptotic order. Because the linear 
operator used to find corrections to the subgrid field is singular — the operator 
necessarily has homogeneous solutions that compose the slow subspace Eq — 
the Fredholm alternative is that the 'right-hand side' of the equation for the 
subgrid fields must lie in the range of the singular operator. This solvability 
condition is enough to determine an extra correction to the evolution. 

Recall from elementary linear algebra that to be in the range of the operator, 
the solvability condition is that the right-hand side must be orthogonal to 
the null space of the adjoint operator. Thus the first task of this section is 
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Figure 7: each element is effectively divided into four subregions by the 
non-locality of the boundary conditions (|9]). To derive the adjoint, label the 
edges of these four subregions as shown. 



to find the adjoint operator of the linear constant diffusion pde ([8j) with its 
insulating boundary conditions (|9]). These non-obvious adjoints have never 
been identified before. Second, we find a basis for the null space. Lastly, 
computer algebra computes the required extra order in the evolution. 

The decoupling of the elements, provided by y = in the boundary condi- 
tions (|9]), simplifies finding the adjoint: we need only consider each element 
in isolation. Thus define the inner product to be the integral over the i, jth el- 
ement: 

rr 

(u, v) uv dx dij = 

To find the adjoint recognise that each element is effectively subdivided into 
four subregions by the coupling of the boundary values to internal values by 
the boundary conditions ([9]), Figure [T] schematically shows these four regions. 
In addition, there exists previously implicit conditions that the subgrid field u 
and its gradient are continuous throughout the element. Then integration by 
parts, or the divergence theorem, transforms the inner product 



Xi- 



uv dx dxj 



(V^u,v) = (u,V^v) 



uv. 



dlj + 



UV- 
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+ 



itxV — uVx dy 



M+ 



+ analogous integrals on B, C and T, 

where specific parts of the boundary integrals are labelled as shown in Figure [7j 
Using superscripts to denote evaluation, continuity requires u*^^ = u*^ and 
vJ^^ — , and the boundary conditions ^ imply u"- — = . Thus 
the inner product 

fYi+1 

V^u,v) = (u, V^v) + 



+ (analogous x integrals of y derivatives) 



= (u,V'v) + 



+ (v^ 



j+1 



u^v^ + u^v^ + uf (v 



M- 



V 



IVH 



) dy 



+ (analogous x integrals of y derivatives). 

For the adjoint, these integrals on the right-hand side must vanish for all 
smooth fields u. Consequently, the null space of the adjoint operator satisfies 
Laplace's equation V^v = with conditions: firstly, that v is zero around 
the edges L, R, B and T of the element; secondly, that v is continuous on the 
interior partitions M and C (but its gradients may be discontinuous there); 
thirdly, that v^-v^"+v^+-v^ = ; and lastly, that v5-vj"+v5^+-vj = . 
Because of these conditions, the null space of the adjoint is spanned by the 
'pyramid' v = (1 — |x — Xi|/H)(1 — |y — Yi|/H) as displayed in Figure Isj 



The solvability condition is then that the integral of the subgrid residuals 
of the PDE with this v over the 1, jth element determines a correction to the 
model evolution. 



It will not escape your notice that the solvability condition integral parallels 
integrals in the Galerkin finite element method. Thus view the Galerkin finite 
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Figure 8: basis 'pyramid' for the null space of the adjoint operator on an 
element: v = (1 - |£,|)(1 - |-n|) = (1 - |x- Xd/H)(l - |y - Yd/H). 



element method as a leading approximation to our systematic slow manifold 
closure of discrete modelling. 

Analogous arguments derive the adjoint for patch dynamics. The version of 
boundary conditions ^ for patches refers to the grid value of the field on 
the right-hand side. With care integrating over a patch, the linear constant 
diffusion operator on the right-hand side of PDE dsj) has adjoint 



3Ei 



9v 

— ds ) 6Dirac(x-Xi,-y-Yj; 



such that V = on 9Ei,j . (20) 



The point source at the central grid point makes up for the diffusive loss of 
material across the edge of the patch to lead to a null vector of the adjoint 
with a logarithmic singularity at the grid point. We have not yet used this 
interesting adjoint for patch dynamics. 
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Return to the discrete Ginzburg— Landau model Computer algebra 
readily computes the subgrid residuals of the Ginzburg-Landau pde ^ to 
the next higher order (the code is detailed elsewhere |34, §2.3]). Taking the 
inner product with the adjoint null vector v, and remembering contributions 
from the inter-element coupling conditions Q, gives the next higher order 
discrete model. Because of the faithful resolution of subgrid structures and 
inter-element interactions, the resulting discrete models are algebraically 
extremely complicated. Thus users may prefer, as they often do now, to 
use models of the nonlinear dynamics of lower order. Then higher order 
discretisations derived via this approach provide local estimates of the local 
error in a lower order simulation as it is computed on the fly. 



4 Generally compute multi-dimensional 
subgrid fields numerically 

Here the subgrid field of the slow manifold is constructed numerically for 
the example Ginzburg-Landau pde ^ in 2D. The approach generalises 
straightforwardly to a wide range of nonlinear pdes in multiple dimensions, 
such as the general reaction-diffusion pde ([T]). 

New complexities arise. Although the spatial structure is obtained numerically, 
the subgrid field must be also parametrised by the grid values Uij, the 
inter-element /patch coupling parameter y, and the nonlinear parameter a. 
Therefore, the construction involves symbolic parameters. The algorithm 
required to develop the holistic model must efficiently solve the corresponding 
mixed numerical and symbolic problem. Only then will the dynamical systems 
methodology be able to be usefully employed in modelling general pdes such 
as ([T]). The focus of this section is on this novel combined algebraic/numerical 
construction of the subgrid fields. 

Numerical construction of the subgrid field introduces errors which are sep- 
arate from the orders of errors in approximating the slow manifold. These 
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errors from the numerical construction of the hohstic discretisation are a ma- 
jor concern. The numerical construction of the subgrid field and its evolution 



has challenging details: ^4.1, How should the subgrid problem be solved? 



^4.2, What subgrid resolutions will accurately reproduce the analytical holistic 



models? ^4.3, What is an efficient implementation? 



This section does not compare the effectiveness of the holistic modelling 
with traditional approaches, that effectiveness has previously been reasonably 
established in ID |32lll2lll8l e.g.]. Instead this section focusses on developing 
the methodology needed to use the approach for more challenging pdes in 
multiple dimensions and in the challenge of supporting the potential efficiencies 
of the gap-tooth scheme of the next section [5j 



4.1 Outline the numerical slow manifold in 2D 

In each element /patch we discretise the microscale subgrid as shown in Fig- 
ure |9] Whereas previously we sought the field Uij(x,'y,t) in the i, jth element, 
here we seek the microscale grid values in the i, jth element /patch. Conse- 
quently, here we index the subgrid by variables k and £ where |k|, \l\ < n. 
For example Figure |9] shows the case n = 4 ; the subgrid is shown solid 
(blue/magenta) for this particular example of a 9 x 9 subgrid. The subgrid 
field extends to ± rH and Yj ± rH to apply either of the 2D non-local 
iCCs Q or ([5]). At the k, £th subgrid point of the i, jth element/patch we 
define the microscale grid value Uij^ic,f (t) for subgrid indices |k|, |£| < n . In this 
section we invoke centre manifold theory to analyse the microscale evolution 
of Uij^i^^f(t) in order to model the macroscale dynamics. 

After discretising the subgrid microscale field in the 1, jth element/patch, 
classic finite differences approximate the spatial derivatives of the subgrid 
field of the Ginzburg-Landau pde (|2]): 



Tony Roberts, January 13, 2013 



4 Generally compute multi-dimensional subgrid fields 



numerically 31 



I 
I 

I 
I 
I 

I . 

1 1- 

I 

-r- 
I 



U4 



1 1 - 
I 

-J-- 



U 1,j 



I 1- 

I 
I 
I 



I 
I 

I 
I 



1+1,3 



-^j - ^ 1,3 - 1 I i+Tj - ^ 



TT 



Figure 9: Example of the n = 4 subgrid on 2D elements (left) and patches 
(right): a subgrid labelled as "n = 4" extends to be 9 x 9 as a microscale 
subgrid extends to in on all sides of a macroscale grid point. 



(21) 



These microscale discretised equations are solved at each of the subgrid points 
inside each of the elements. The overlapping elements are coupled with iccs 
corresponding to Q, namely 



Ui,j,±n,« = yui±i,j,o,i! + (1 - y)u-y,o,f , 1^1 < n , 

_Ui,j,k,±n = yUi,j±i,k,0 + (1 - y)Ug,k,0 , M < TL , 

or in the case of disjoint patches by iCCs corresponding to ([s]), namely 

'^,],±n,i = £t^iy]£f^'^[y]T^i,]fifi > 1^1 < , 
uy,k,±n = ^^[^^'^(y)^:^ '■(y)ui,j,o,o , |k| < n . 



(22) 



(23) 



Equations (21)-(22)/(23) form a system on a microscale lattice. The same 



centre manifold theorems apply to the system (21)-(22)/(23) to assure us 
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of the existence, relevance and construction of a slow manifold, macroscale 
discrete model of the lattice dynamics. Indeed the application of centre 
manifold theory is more straightforward as the microscale lattice system is 
finite dimensional, in contrast to the 'infinite' dimensionality of a microscale 
PDE. The macroscale slow manifold to construct is that the subgrid field 
■Ud,j,k,« = Uy^i,,f(U,y, a) where, defining the grid value Uy = Uy^co , the 
macroscale grid values U evolve according to the system U^j = gij(U,y, oc). 

An iteration scheme finds the microscale subgrid field and the macroscale slow 
evolution. We leave the case of patches until Section [5] and here focus upon the 
case of overlapping elements. The initial approximation is that of a constant 
field in each element: ik^jxi ~ such that U^j = g^j pa 0. Given any 
current approximation, Uij and gij, we seek an improved approximation to 
the field, uij +^1^, and evolution, g^j + g^, where Uy and g^ are corrections 
to be found in each iteration. At each iteration, and in the case of overlapping 
elements, the following linear equations driven by the current residuals are 
solved for the corrections Ug and g( j 

Kj,±n,f - W,j,o,f = ReE|22], ul]x±n - W,j,k,o = Re^, Ug = • (24) 

The iteration repeats until all residuals are zero to a specified order of error. 
This iteration scheme follows that for the analytic construction of the subgrid 
field and is documented in full detail in a separate technical report §4]. 
Centre manifold theory then assures us that the resultant macroscale discrete 
model Uij — gij(U,y, oc] is accurate to the same order of error. 

For example, for the coarsest possible subgrid field, the case n — 2 when the 
microscale subgrid within each element has half the spacing of the macroscale 
grid, the lattice dynamics on the microscale subgrid gives the low order, 
macroscale model 
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^6^Ug + ay - m,8'U,,) + 0{y' + oc') , (25) 



Compare this model with the analytic macroscale discrete model (15): the 
terms in the first line are identical; the same higher order terms appear in 
the second line but the coefficients differ by 25%. This correspondence is 
promising for such a coarse microscale subgrid discretisation. 

Recall that this approach models the dynamics on a microscale lattice. Thus 
this approach transforms microscale lattice dynamics onto a macroscale lattice. 
Such transformation across scales of lattices may be repeated across an entire 
hierarchy of lattices to form dynamics on a multigrid [17]. The consistency 
theorems of section [6] provide additional suuport for the modelling across a 
hierarchy of lattices. 



4.2 Low resolution subgrids are accurate in 2D 

How do macroscale models constructed via a numerical microscale, such 



as (25), compare with analytic macroscale models? We compare in two ways: 
one via the convergence of the coefficients; and the other by the accuracy of 
the predicted bifurcation diagrams. It appears that the microscale subgrid 
need not be of high resolution. 



First look at the coefficients of the 0(y^ + cc^) models such as (25). Recall 
that the number of microscale subgrid points, from one macroscale grid point 
to the next, in each dimension, is n. The coefficients linear in the coupling 
parameter y and nonlinearity a are exact for all n > 2 , only the higher 
order coefficients vary with subgrid resolution. Table [T] tabulates coefficients 



in these nonlinear terms, those of 0(y^ + a^) in models such as (25), for 
some values of n. Evidently the coefficients converge to the continuum pde 
values with error 0(l/n^). We expect such quadratic convergence from the 
quadratic approximation in ([2T]) of the subgrid scale PDE dynamics. 



Similarly, we compare numerically obtained coefficients for the O (y'^ + oc^) 
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Table 1: coefficients in the (9(y^ + a^) models, such as (25), evidently converge 
to the correct analytic coefficients, in (15) and labelled oo in the table, with 
errors /n^) as the resolution of the microscale grid improves. 

...u„.-j -.2c4ir /111 ,„.c2ii3 i2 c2i i 



2 


5x5 


1 

16 


1 

16 


3 
16 


4 


9x9 


5 

64 


5 

64 


15 

64 


8 


17 X 17 


21 

256 


21 

256 


63 
256 



OO 



_1_ 

" 12 



_1_ 
12 



Table 2: maximum errors in the coefficients of the 0(y'* + Oi^) model when 
approximated numerically at three different subgrid resolutions. The decrease 
by at least a factor of four, upon doubling n, indicates quadratic convergence, 
n y^/H^ ya, y^/H^ y^oc ya^H^ 
2 0.021 0062 00033 014 OOOTe 
4 0.0052 0.016 0.00086 0.040 0.000098 
8 0.0013 0.0039 0.00022 0.010 0.000006] 



model. Because of the complexity of the model we make a limited comparison: 
for each order in a and y , Table [2] reports the largest error in the numerically 
obtained coefficient for three different subgrid scale resolutions. Evidently, 
these maximum errors decrease like 1/n^ to confirm the accuracy of the 
numerical description of the subgrid scale dynamics. 

Second, we turn to the bifurcation diagram to see the sort of errors incurred 



in using the approximate models. Figure 10 shows the bifurcation diagrams 
for the (9(y^, cc^) holistic model of the Ginzburg-Landau system for four 
subgrid resolutions. Here the equilibria shown in green are not the accurate 
solution of the Ginzburg-Landau system, but rather the equilibria of the 



analytic 0(y^, oc^) holistic model in 2D (obtained from (15) by omitting the 



y term). Observe that with a subgrid resolution of just 4x4 intervals the 
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(a) 2D-RGL 0{y^,a^), 2x2 interval subgrid 
8 




(c) 2D-RGL 0(/,a^), 6x6 interval subgrid 
8 




(b) 2D-RGL 0{y^,<x), 4x4 interval subgrid 
8 




(d) 2D-RGL 0(Y^,a^), 8x8 interval subgrid 
8 




Figure 10: Bifurcation diagrams of the 0(y^, 0(}] holistic models of the 2D 
Ginzburg-Landau system with 8x8 macroscale elements on [0, 7t] x [0, tt] for 
subgrids of (a) 5 x 5, (b) 9 x 9, (c) 13x13 and (d) 17 x 17. The bifurcation 
diagram for the analytically constructed model is shown in green. 
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bifurcation diagram for the numerically constructed 0(y^, (X^) holistic model 
is almost indiscernible from the analytic model over nonlinearity < a < 20 . 
Higher subgrid resolutions are indiscernible to even larger nonlinearity a. 

As a last comparison of bifurcation diagrams, Figure [6] shows one example 
confirming that by computing to higher order in the inter-element interactions, 
and resolving the subgrid scale structures numerically, the predictions of the 
numerically derived models do improve significantly. 

Numerically resolving the microscale subgrid structures does generate usefully 
accurate, slow manifold, macroscale, discrete models. 



4.3 An efficient computer algebra approacli is crucial 

The difficulty associated with the numerical construction of the subgrid field 
is the mixed numerical and symbolic nature of the equations involved in 
the iteration scheme. The size of the system of equations increases as the 
subgrid resolution improves, and the complexity of the symbolic nonlinear 
residuals increases quickly as higher order holistic models are constructed. 
However, these comments apply to, and this subsection only addresses, the pre- 
simulation preprocessing step of constructing the slow manifold discretisations: 
the computation time of the derived macroscale closure on the macroscale 
domain is as normal for the method of lines, and is independent of the 
considerations here. 

Computer algebra packages such as Reduce \18\ or Mathematica [56] have 



general routines that solve systems of equations such as (24). However, 
these solve routines are inefficient for the many equations and complicated 
expressions here. Even a low resolution, n = 2 subgrid, took many minutes 
using the built in solve routines of both Reduce and Mathematica. Instead 
we develop an approach that is practical for implementation with a large 
number of complicated symbolic terms. Nonetheless, different computer 
algebra packages execute this preprocessing step in a time that differs by an 
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Table 3: Reduce and Mathematica computational times for numerical construc- 
tion of 0(y'^, (X^) holistic models of the one dimensional Ginzburg-Landau 
equation for various subgrid scale resolutions, n. 



n 


Reduce 


Mathematica 


2 


1.1 s 


70.2 s 


4 


3.1 s 


215.4s 


8 


8.3 s 


367.6 s 


16 


23.7s 


517.7s 



order of magnitude: users need to know that so far we find Reduce [T8] to be 
the most efficient. 



Transform to constant coefficient Recall that at each step of the iter- 
ation scheme we solve a problem for updates to the subgrid field u(j and 



its evolution g[y Multiply the first (field) equation in (24) by r H and 



replace g(j by Q' — r^H^g(j . Then the left-hand side of the new form of 
the equations has numerical constant coefficients; algebraic expressions only 
occur in the right-hand side. 

Further, the left-hand side of the new equations remain the same for every 
iteration. Consequently, the first iteration constructs an LU factorisation of 



the left-hand side, which is then used to solve equation (24) for updates in 
every iteration |3H §4]. The LU factorisation is performed once and requires 
approximately ^N^ operations |10| e.g.]. Here the number of equations for the 
subgrid structure are N = (2n + 1 )2 + 1 ; for example, N = 290 for the n = 8 
microscale subgrid. At each step of the iteration scheme the LU factorisation 
algorithm operates on the symbolic residual vector. We suspect 2D and 3D 
problems would be solved more efficiently through sparse methods such as 
iterative multigrid |36[ [2] or incomplete LU factorisation and Krylov subspace 
methods [231 [51]. Such alternatives remain for later exploration. 
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Reduce was much faster Computational experiments found that the 
computer algebra package Reduce was at least an order of magnitude faster 
than Mathematica. Table [3] lists the computational time for the Reduce 
and the Mathematica implementation for constructing 0(y'^,<x^) holistic 
models of the one dimensional Ginzburg-Landau equation with subgrid 
resolutions of 2, 4, 8 and 1 6 subgrid intervals. These times were observed on a 
Pentium III, 750 MHz processor, with 256 Mb RAM, running Reduce 3.7, under 
Windows XP. Table |3] shows the Reduce implementation was 20-70 times 
faster than the Mathematica implementation (even with the repeated help of 
the Mathematica news group). Thus we use the free package Reduce [18] . 

5 The slow manifold of emergent patch 
dynamics 

Now return to the gap-tooth dynamics on patches such as the simulations 
of the Ginzburg-Landau pde ^ shown in Figures [l}|3j The challenge is to 
deduce by analysis some of the important properties of such 'equation free' 
simulation so we understand its performance in applications. 

This section shows how the emergent patch dynamics may be viewed as a 
slow manifold closure that is also consistent with the underlying microscale 
model. Hence, for example, patch dynamics could explore the interesting 
domain competition inherent in the Ginzburg-Landau pde. 

Subgrid microscale modes form very rapid transients Figures [l}]3] 
sample the evolution from a non-smooth initial condition of the Ginzburg- 
Landau PDE ([2]). The transients are the rapid diffusive smoothing of the 
initial jagged microscale structures within each spatial patch. Thereafter we 
would see the slow macroscale evolution of domain competition, underpinned 
by smooth microscale structures within each patch. Thus in the long term 
we see relatively slow evolution of smooth structures in each patch. 
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100 
decay rate 



200 



Figure 1 1 : histogram of the number of eigenvalues of linear diffusion — oc — O 
in the Ginzburg-Landau pde ^ — showing clear time scale separation between 
the 1 6 global modes with small eigenvalues, and the multitude of microscale 
modes decaying at rates > 40. Here the gap ratio r = 1/4 and macroscale 
spacing H = 2. 
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Eigenanalysis of linear diffusion — <x = in the Ginzburg-Landau pde ^ — 
confirms this separation of modes into rapidly decaying subgrid modes and 
long lasting macroscale modes. For example, consider doubly periodic diffusion 
on 4 X 4 patches on a macroscale grid with spacing H = 2 with each patch 
composed of a 5 x 5 microscale grid. Couple patches with the iCCs ([s]). Then 
numerical differentiation of the code that simulates Figures [l]-[3] generates 



the matrix of the diffusion dynamics on these patches. Figure 11 plots a 
histogram of the eigenvalues of this matrix. Most eigenvalues are large and 
negative corresponding to the rapid diffusive decay of subpatch structures. 
However, 16 eigenvalues are near zero, one for each of the 16 patches. These 
near zero eigenvalues characterise long lasting, emergent, global modes. 



Figure 12 confirms the smooth global modes. Not discernible in the eigenvalues 



of Figure [TT] is that the small eigenvalues form four main groups. Figure 
plots a representative eigenmode from each of these four groups: from the 
near sinx mode of A = —0.471 to the macroscopic zig-zag but microscopic 
smooth mode of A = —2.05 . This linear analysis confirms the basis for 
theoretical support of nonlinear patch dynamics for general reaction-diffusion 

PDES 



Centre manifold theory supports macroscale closure Section |2] in- 
troduces y to parametrise the coupling between patches: y = 1 corresponds 
to the coupling used in 'equation free' simulation, but when y = the eigen- 



values are perturbed so that the small eigenvalues, such as those in Figure 11 



become precisely zero. Then Corollaries [T] and |2] establish the existence of 
emergent patch dynamics as a slow manifold at finite coupling y and finite 
nonlinearity a for general reaction-diffusion pdes ([T]). 



Patch dynamics are consistent with the microscale Section |4] de- 
scribes how to construct slow manifold discretisations of general reaction- 
diffusion PDES ([T]). For example, consider the Ginzburg-Landau pde ^ when 
discretised inside patches by a (2rL+ 1 ) x (2n+ 1 ) microscale lattice. Using the 
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> =-0.488 > =-0.986 




Figure 12: some global eigenmodes labelled by their eigenvalue. Here the 
ratio r = 1/4 and macroscale spacing H = 2. For comparison, a second 
order, finite difference scheme with corresponding grid spacing H = 2 would 
have eigenvalues A = — 1 /2, — 1 , — 3/2, — 2: the shown macroscale eigenvalues 
correspond well with these finite difference ones. 
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patch iccs (|5| computer algebra computes the slow manifold discretisation is 



y 



(1-&)6Xi + ^(l-S) 



12H2^ n^/ ^0 90H2 



4n2 



(26) 



for operator 6 defined by ( 16 ) . The nonstandard discretisation of the cubic 
reaction, appearing via the ccy term above, is obtained from a systematic 
approximation of the subgrid, microscale, out of equilibrium, structures: its 
coefficient varies with the {2n + 1 ) x (2n + 1 ) microscale lattice as 



0.0758188 as TL ^ oo , and Cr 



751 - 
6(1651 + 191/n2 



exactly for n = 2, 3, 4 and to four significant digits for n = 5 : 8 . Evaluated at 
full coupling y = 1 corresponding to the 'equation free' gap-tooth simulation, 
the discretisation (26) is consistent with the microscale lattice discretisation 
of the Ginzburg-Landau pde (|2]). Replacing the differences in (26) by their 
expansion in derivatives, the equivalent pde to (26) is 



9U 



2^2 



12n2 V9x4 ^ 



+ (x(U - U^) + cxHh^eCnU 



4^,4 



360n4 
9U\2 



a^u a^u 

ax^ dy^ 

auy 



(27) 



This PDE is not the Ginzburg-Landau pde ^ because the underlying dynam- 
ics are those of the microscale discretisation: instead this pde is equivalent to 
that of the microscale discretisation on the fine lattice which has spacing rH/4. 
Such consistency of macroscale patch dynamics with the fine scale dynamics 
is established for a range of pde operators in Theorem [6] of the next Section |6j 

Lastly, these patch dynamics connect to classic finite difference discretisa- 
tions. Fix the macroscale grid spacing H, but let the patch size and the 
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microscale grid spacing become small via the ratio r — > . Then the slow man- 



ifold model (26) reduces to the classic finite difference approximation of the 
Ginzburg-Landau pde ([2]): the equivalent pde (27) is the Ginzburg-Landau 
PDE ([2]) as patch size r — > 0. Using classic interpolation to provide coupling 
conditions ^ for tiny patches is equivalent to classic finite differences of 
the microscale dynamics. Thus when we apply the 'equation free' gap-tooth 
method on simulators for which we do not know the microscale equations, we 
will nonetheless obtain a macroscale simulation which is consistent with the 
unknown microscale equations. 

Furthermore, this consistency of the macroscale simulation holds no matter 
how small the patches, here parametrised by the ratio r. The constraints 
on the macroscale grid will be those familiar to anyone using classic finite 
difference or finite element approximations: namely, the macrogrid step H 
must be small enough to resolve the macroscale variations. But the patch size 
can be vanishingly small making the gap-tooth method extremely efficient for 
those systems with extremely wide separation of scales. 



6 Non-local coupling conditions enforce 
consistency 

Recall that the constructed holistic models of the Ginzburg-Landau pde ^ 
are consistent with the pde as the grid size H ^ 0, see equations (18) in 
Section [3| Now we prove that general consistency follows from the specific 
choice of nonlocal inter-element/patch coupling conditions ^ and ([s]). 

We start with a similar theorem to one previously proved for the consistency 
of holistic discretisation in one space dimension [1^ . The critical innovation 
here is in the proof: previous proofs were constructive whereas here it is not. 
Avoiding a constructive proof is essential here as we do not know analytic 
forms for the slow manifold subgrid fields in multiple dimensions. In this new 
proof: a new lemma establishes consistency for nonlinear reaction-diffusion 
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PDEs in ID; and an immediate new corollary then proves consistency of the 
2D holistic discretisation of a wide range of 2D nonlinear, reaction-diffusion 

PDEs. 

Theorem 3 (ID linear consistency) Consider the pde 9tU = £u for 

some local, isotropic, homogeneous, linear operator C Model the dynamics 
on overlapping elements of an equi-spaced grid Xt = iH . Let it^(x, t) denote 
the subgrid field in the ith element satisfying the pde QtXti — Cxii on the 
interval (Xi_i,Xi+i) with the moderated inter-element coupling conditions 

Ui(X^±i,t) =yit^±i(X,±i,t) + (1 -y)ui(X,,t). (28) 

When inter-element interactions are truncated to residuals O(y^) the grid 
values Ui(t) = Ui(Xi, t), at full coupling y = 1 , evolve consistently with the 
PDE 9tU = Cu to an order in grid size H that increases with p. 

Proof: We proceed with some classic operator algebra PB] e.g.]. The 
principle obstacle is to transform subgrid spatial differences, indicated by 
subscript x, into macroscale grid differences, indicated without a subscript. 
Begin with the pde on the ith element: 9tUi — Cui . Because the operator C 
is isotropic and homogeneous it may be formally expanded in even centred 
differences as 



k=0 



for some coefficients lik and corresponding function £. For example, in 
application to reaction diffusion PDEs, we expand the diffusion operator a^ 



9x = n^i) 



^sinh-^(l5,] 



1 



= — \8^- , ±56 _ _L58 ... 1 



^Such operator expansions and our formal operator manipulations appear to be little 
known these days, but they are well established [38l|T9l e.g.]. The manipulations are valid 
for infinitely differentiable functions as appropriate to the diffusion-like pdes considered 
here, but not applicable to systems generating shocks or other singularities that are not the 
subject of this theorem. Modern analysis typically prefers to invoke Taylor's Remainder 
Theorem which avoids requiring infinite differentiability, but the aim here is to prove 
consistency to arbitrarily high order. 
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which has the inverse 

62 = r^(92) ^4sinh^(H3,/2) 



In general, provided the leading coefficient ^2 7^ , the pde dtUi = [^o+^(5^)]ui 
is equivalently written l^^[dt — ^ojiJ-i = 5^ui where is the inverse of 
function £ (since ^ , I is invertible for small enough argument). 

Evaluate this last form of the pde at x = so that the Ut on the left-hand 
side becomes simply the grid value Ui and by the coupling conditions (28]^ 
the centred spatial difference on the right-hand side becomes the centred grid 
difference yS^Ui . This evaluation then gives the evolution to be £^V9t — 
£o)Ui = yS^Ui on the macroscale grid. 

Now reverting the inverse function, this grid evolution is equivalent to 

oo 

9tU, = [£o + ^yS')] U, = Y_ W'S^'^U, . (29) 

k=0 

For example, for the diffusion equation 

9tU, = — [2sinh-i (lVyS)]'Ut 

:5 + 



1 
1 



^ M 90 560 



Thus a truncation of (29) to errors O(y^) results in a discrete model with 
stencil width of 2p — 1 . But specifically relevant to the theorem is the 
equivalent differential equation of this discrete model evaluated at full coupling: 
for smooth macroscale U^, the error in approximating C by the truncated 
version of (29) (at full coupling y = 1) is dominated by the leading neglected 



term, namely ^ipS^^. As the element size H ^ 0, this error is (^(^ipH^P) for 



^If the leading coefficient in the expansion of £ is l2n , because the lower order 
coefficients are zero (or asymptotically small as in the Kuramoto-Sivashinsky pde), then 



more coupling conditions like (281 couple with the next nearer neighbouring elements. 
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smooth fields. For example, for the diffusion operator 9^, the coefficients 
£2k = 0(^ /H^) and so the discrete model is consistent with the diffusion pde 
to error 0(H^^^^) as grid size H ^ . ^ 



The above proof is so slick that one is tempted to immediately apply the 
identical arguments to nonlinear pdes. However, the stumbling block is 
that in the reversion of the operator £ we need time derivatives and spatial 
differences to commute with nonlinearity. Such commutation is not exact, 
only approximate, and leads to the following modification of the proof of 
consistency to nonlinear pdes. 

Lemma 4 (nonlinear consistency) Consider the nonlinear, reaction dif- 
fusion PDE 9tU — £u+g(u) for some local, isotropic, homogeneous, conserva- 
tive, second order, linear operator C and some smooth nonlinear reaction g(u). 
Model the dynamics on overlapping elements of an equi-spaced grid Xt = IH . 
Let iti^(x, t) denote the subgrid field in the ith element satisfying the PDE 
QtXii — £ui + g[vii) on the interval (Xt_i,Xt+i) with the moderated inter- 



element coupling (28). When inter-element interactions are truncated to 
some nonlinear order in y, the grid values Ui,(t) = Ui,(Xi,t), at full coupling 
y — } , evolve consistently with the PDE 9tU = £u+ g(u). 

Proof: Slightly different to before, the operator C may be formally expanded 
in even centred differences as 

1 °° 1 

k=l 

for some coefficients Izk and corresponding function £. But here: firstly, since 
C is conservative, the coefficient = as here the reaction g(u) absorbs any 
non-zero Iq; and secondly, since the operator £ is second order, = C'(l) 
as H — > ; the main H dependence is explicit in the above expression. 
Noting the inverse £-i(X) = l-X + 0{X^) and the PDE in each element is 
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H^9tUi = £(5^)ui + H^g(ui), we observe 

{H%]ui = ^H%Ui + O (H^) 



1 



= 6^ut+— g(u,) + 0(H4), 
''I 

when the fields u and Ui and reaction g are smooth enough so that their 
(2k)th order differences are 0(H^'^) and time derivatives are 0(^) as H ^ . 
Evaluate the above equation at the macroscale grid points x — Xi, using the 
coupling condition ( 28 ) , to determine 



{H^d,)Ui = yS^Ui + ^g(Ui) + ^(H^) . 
Use the approximation to the operator to obtain 



yS^Ui + ^glUO 



+ C(H^) from g 



£(y6^)Ut + H^g(U0 + O(H4). 



Thus, dividing by 



(30) 



which, for any truncation of O(y^) or higher and then evaluated at full 
coupling y = 1 , is consistent as H ^ with the nonlinear reaction-diffusion 
PDE itt = £u + g(u). 4 
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We expect that more careful treatment of the O (H^) error terms in the above 
will show the error is of higher order in the grid size H, as observed in all 
examples, but we do not attempt to do so here. 

In this section the subgrid microscale operator C need not be differential. 
For example, C could be a microscopic lattice operator as in the numerical 
construction of the previous Section |4j for example, with n microscopic grid 
points for every macroscopic grid point we have the microscale lattice diffusion 

6^;^^,^ = 4sinh^ [\hd^] = 4sinh^ 

The proof still holds. Further, the proof still holds if the subgrid field is 
not just a scalar. This last observation empowers the following corollary on 
general consistency in two dimensions that Section |3] observed specifically for 
the Ginzburg-Landau pde. 

Corollary 5 (2D consistency) Consider a reaction- diffusion pde 9tU = 
V^u+g(u) (such as the Ginzburg-Landau equation ^) modelled on overlap- 
ping elements with subgrid fields u^[x,y,t) coupled by conditions When 
the interactions are truncated to some nonlinear order in y the grid values 
Uy(t) = Uij(X^, Yj,t), at full coupling y = 1 , evolve consistently with the 

PDE. 

Proof: Apply the previous Lemma |4] twice. First, treat coordinate y as 
a parameter so that the reaction is 0(u) = g(u) + d^u and ii^l) = Q^- 
Then by Lemma |4] the semi-discrete 'grid' values XLij(Xi,'y,t), for discrete I 
and parametrised by continuous y, evolve consistently with the reaction- 
diffusion PDE. Second, treat index i as a parameter and consider the discrete 
modelling in coordinate y so that now the reaction g also involves operators 
acting on the x-grid and = 9^. Then by Lemma |4] the 2D grid values 
Xli j = Ui^j [Xi, Yj , t) evolve consistently with the semi-discrete system generated 
in the first step, which in turn is consistent with the diffusion PDE. ^ 

The generalisation to higher spatial dimensions appears immediate. 
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The following adaptation to the dynamics on patches is weaker, but the 
theorem is a starting result to support the efficient equation-free modelling 
on patches [52l e.g.]. As centred differences generally scale like 6 oc H 
as macrogrid spacing H ^ 0, the theorem confirms simulation errors will 
be 0{H'^'^^^) where in the exponent I characterises the system being simulated 
(1 = 2 for a diffusion operator). 

Theorem 6 (isotropic patch consistency) Consider the pde 9tU = Cu 

for some local, isotropic, homogeneous, linear operator C in 2D. Model the 
dynamics in patches of size ft = rH centred on an equi-spaced grid Xi = IH 
and Yj = jH . Let uyfxjiijt) denote the suhgrid field in the (i, patch 
satisfying the pde StU^j = >Cuij on the patch with the moderated interpatch 
coupling conditions ([S]). When interpatch interactions are truncated to resid- 
uals 0{y^^ the grid values Uy(t) = Uij(Xi,, Yj,t), at full coupling y = 1 , 
evolve consistently with the pde 9tU = £u to errors 0{b^ + 6^^) . 

Proof: First let's establish how patch sized differences relate to grid value 
differences. Second, this relates the pde with the evolution of grid values. 
Define the patch sized centred difference by/\x[x) = u[x + 1t./2) — u(x — fi/I). 
Then evaluated at the grid point (X^, Yj) the second centred difference 

K^j] (X„Yi) = [{^x - 2 + f rKj](X„Y,) 

= {8l[y)-l + er' (y ) }Ui,j by iccs Q , and then by Q 
= { [1 + y(Mi + 2 + [1 + y(-Mi + Ui,j 

= { [1 + y(^,5, + \b\]Y- 2 + [1 + y(-^,6t + \b\)Y] Uy 

+ (9(y^5^^) upon truncating to errors 0(y'') 
= { [1 + ^itSt + 16?]' - 2 + [1 - ^tSt + 16?]'} Uy 

+ (9(6?^) upon evaluating at y = 1 

= {f];-2 + £:r}Uy + 0(6?) 

= S?U,, + (9(6?) 
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Similarly, the analogous relation holds in the other spatial direction: 

Now relate the PDE and the patch dynamics. Because the linear operator C 
is isotropic we may formally write it in terms of the isotropic operator as £ = 
lQ + l[8^ + 5y). Then inverting £ the system of pdes on patches, 9tUi,j = /^Udj , 
becomes f^^Qt — ^ojiJ-y = i^i + ^y)'^,j ■ Evaluating at the grid point (X^, Yj) 
and using the above results leads to (at-^ojUy = (S^+Sf)Uy+C(6f +6^'') . 
Lastly, reverting I gives dtU^^j = [to + + Sf)]Uy + 0{bf^ + 6^^) = 

CUi^j + 0(5^^ + 5^^). That is, the grid values evolve consistently with the 
system simulated within the patches. 4 



7 Conclusion 

We explored novel macroscale discretisation of reaction-diffusion dynamics in 
two spatial dimensions. This work generalises considerable earlier work on 
modelling dynamics in one dimension. The specific coupling conditions ^ 
and ([5]) have important theoretical and practical consequences for inter- 
element and interpatch dynamics 

Section [2] discussed how these coupling conditions ensure that centre manifold 
theory applies to prove the existence, emergence and approximation of the 
slow manifold that is the macroscale discretisation of general reaction-diffusion 
equations in two dimensions. Further, Section [6] proved that the resulting 
discrete models will also be consistent, as the macroscale grid size H ^ , 
with the continuum or microscale dynamics. Thus the holistic discretisations 
generated in this novel approach have the dual justification of both consistency 
for small H and the existence, by centre manifold theory, of an exact slow 
manifold at finite H. 
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This strong theoretical support appears to be straightforwardly generalisable 
to reaction-diffusion dynamics in three or more dimensions. The support also 
appears to be straightforwardly generalisable to higher order pdes, just as 
the theory supports the discrete modelling of one dimensional higher order 
PDES such as the Kuramoto-Sivashinsky pde [301 132]. 

Sections [3] and |4| using the example of the real Ginzburg-Landau pde (|2| , 
explored technical issues necessary to apply the approach here and to more 
general PDEs. In contrast to one spatial dimension, a purely algebraic approach 
can only carried out to a low order of accuracy. Consequently we introduced 
and explored an approach where the microscale subgrid dynamics are described 
numerically, but with algebraic expressions for coefficients so that we construct 
an algebraic model for the macroscale discretisation. This work provides a 
powerful approach and theory for the sound and accurate closure of macroscale 
simulations. 

Section[5]showed that the adaptation of this approach to the gap-tooth method 
of Kevrekidis et al. f2M E2] is sound, even for very small patches of simulators. 
Thus very efficient simulations are possible. However, further research needs 
to extend this argument to systems where the microscale dynamics are not 
the spatially smooth dynamics of reaction-diffusion equations ([T]) . 

Acknowledgement The Australian Research Council Discovery Project 
grants DP0774311 and DP0988738 helped support this research. We thank 
Yannis Kevrekidis for his encouragement and many stimulating discussions. 



References 

[1] Vemuri Balakotaiah and Hsueh-Chia Chang. Hyperbolic homogenized 
models for thermal and solutal dispersion. SI AM J. Appl. Math., 63:1231- 
1258, 2003. |doi:10.1137/S0036139901368863[ 



Tony Roberts, January 13, 2013 



References 



52 



[2] W. L. Brigg, V. E. Hanson, and S. F. McCormick. A multigrid tutorial. 
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 
2000. 

[3] J. Carr and R. G. Muncaster. The application of centre manifold theory 
to amplitude expansions. II. Infinite dimensional problems. J. Diff. Eqns, 
50:280-288, 1983. 

[4] Zhiming Chen and Thomas Y. Hou. A mixed multiscale finite element 
method for elliptic problems with oscillating coefficients. Math. Comp., 
72:541-576, 2002. 

[5] C. Chi cone and Y. Latushkin. Center manifolds for infinite dimensional 
nonautonomous differential equations. J. Differential Equations, 141:356- 



399, 1997. http : / /www . ingentaconnect . com/content/ap/ de/ 1997/ 



[00000141/00000002/art03343: 
[6] Alexandre J. Chorin and Panagiotis Stinis. Problem reduction, renormal- 



ization, and memory. Technical report, |http://arXiv.org/abs/math 
|NA/0503612|, 2005. 

[7] E. J. Doedel, R. C. Paffenroth, A. R. Champneys, T. F. Fairgrieve, Yu. A. 
Kuznetsov, B. Sandstede, and X. Wang. Auto 2000: Continuation and 
bifurcation software for ordinary differential equations (with HomCont). 
Technical report, Caltech, 2001. 

[8] J. Dolbow, M. A. Khaleel, and J. Mitchell. Multiscale mathematics 
initiative: A roadmap. Report from the 3rd DoE workshop on multiscale 



mathematics. Technical report. Department of Energy, USA, http 
[//www .sc. doe ■ gov/ascr/mics/amr , December 2004. 



[9] Weinan E, Bjorn Engquist, Xiantao Li, Weiqing Ren, and Eric Vanden- 
Eijnden. The heterogeneous multiscale method: A review. Tech- 



nical report, |http : // www . math . princeton . edu/mult iscale/review 

[pdf , 2004. 



Tony Roberts, January 13, 2013 



References 



53 



[10] Shin-Ichiro Ei, Kazuyuki Fujii, and Teiji Kunihiro. Renormahzation- 
group method for reduction of evolution equations: invariant manifolds 
and envelopes. Annals of Physics, 280:236-298, 2000. 

[11] B. Ermentrout. XPPAUT 5.0 - the differential equations tool. Technical 
report, ^http : / /www . math . pitt . edu/~bard/bardware/ xpp_doc . pdf j , 
2001. 

[12] Martin J. Gander and Andrew M. Stuart. Space-time continuous analysis 
of waveform relaxation for the heat equation. SIAM Journal on Scientific 
Computing, 19 (6): 20 14-2031, 1998. 

[13] C. W. Gear, Ju Li, and I. G. Kevrekidis. The gap-tooth 
method in particle simulations. Phys. Lett. A, 316:190-195, 2003. 
|doi:10.10l67j.physleta.2003.07.004 



[14] J. D. Gibbon. Weak and strong turbulence in the complex Ginzburg- 
Landau equation. In G. R. Sell, C. Foais, and R. Temam, editors. 
Turbulence in fluid flows-A dynamical systems approach, volume 55 of 
The IMA volumes in mathematics and its applications, pages 33-48. 
Springer- Verlag, 1993. 

[15] H. S. Greenside. Spatiotemporal chaos in large systems: the scaling 
of complexity with size. Technical report, http : / / arXiv . org/ abs/| 
|jchao-dyn/9612004, December 1996. 

[16] H. S. Greenside and W. M. Coughran. Nonlinear pattern-formation near 
the onset of Rayleigh-Benard convection. Phys. Rev. A, 30:398-428, 
1984. 

[17] Bjorn Gustafsson and Jacqueline Mossino. Non-periodic explicit homog- 
enization and reduction of dimension: the linear case. IMA Journal of 



Applied Mathematics, 68:269-298, 2003. doi:10.1093/imamat/68.3.269 



[18] A. C. Hearn. Reduce. Technical report, [http : //vi-wvi . reduce- algebra 
co^, 2004. 



Tony Roberts, January 13, 2013 



References 



54 



[19] Francis B. Hildebrand. An introduction to numerical analysis. Dover, 
1987. 

[20] Thomas Y. Hou and Xiao-Hui Wu. A multiscale finite element method for 
eUiptic problems in composite materials and porous media. J. Comput. 
Phys, 134:169-189, 1997. 

[21] J. M. Hyman, B. Nicolaenko, and S. Zaleski. Order and complexity in the 
Kuramoto-Sivashinsky model of weakly turbulent interfaces. Physica D, 
23:265-292, 1986. 

[22] M. S. Jolly, I. G. Kevrekidis, and E. S. Titi. Approximate inertial mani- 
folds for the Kuramoto-Sivashinsky equation: analysis and computations. 



Physica D, 44:38-60, 1990. doi:10.1016/0167-2789(90)90046-R. 



[23] C. T. Kelley. Iterative methods for linear and nonlinear equations, 
volume 16 of Frontiers in applied mathematics. Society for Industrial 
and Apphed Mathematics SIAM, 1995. 

[24] I. G. Kevrekidis, G. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Run- 
borg, and K. Theodoropoulos. Equation-free, coarse-grained multiscale 
computation: enabling microscopic simulators to perform system level 
tasks. Comm. Math. Sciences, 1:715-762, 2003. 

[25] loannis G. Kevrekidis and Giovanni Samaey. Equation-free multiscale 
computation: Algorithms and applications. Annu. Rev. Phys. Chem., 
60:321—44, 2009. 

[26] Y. A. Kuznetsov. Elements of applied bifurcation theory, volume 112 of 
Applied Mathematical Sciences. Springer- Verlag, 1995. 

[27] C. D. Levermore and M. Oliver. The complex Ginzburg-Landau equation 
as a model problem. In P. Deift, C. D. Levermore, and C. E. Wayne, ed- 
itors. Dynamical systems and probabilistic methods in partial differential 
equations, volume 35 of Lectures in Applied Mathematics, pages 141-190. 
American Mathematical Society, 1996. 



Tony Roberts, January 13, 2013 



References 



55 



[28] Yoram Louzoun, Sorin Solomon, Henri Atlan, and Irun R. Cohen. Mod- 
eling complexity in biology. Physica A, 297:242-252, 2001. 

[29] D. H. Lowndes, D. 0. Welch, 1. McNulty, A. P. Alivisatos, R. S. Averback, 
A. J. Nozik, M. Alper, T. A. Michalske, J. A. Eastman, and T. P. 
Russell. Nanoscale science, engineering and technology research directions. 
Technical report, U.S. Department of Energy, 1999. 

[30] T. Mackenzie and A. J. Roberts. Holistic finite differences accu- 
rately model the dynamics of the Kuramoto-Sivashinsky equation. 



ANZIAM J., 42(E) :C918-C935, 2000. http : //anziamj . austms . org 
[au7V427CTAC99/Macki 

[31] T. MacKenzie and A. J. Roberts. Holistic discretisation of shear 
dispersion in a two-dimensional channel. In K. Burrage and 
Roger B. Sidje, editors, Proc. of 10th Computational Techniques 
and Applications Conference CTAC-2001, volume 44, pages C512- 
C530, March 2003. http : / / j ournal . austms . org . au/ o j s/ index . php/| 
pZIAMJ/article/view/694, 

[32] T. MacKenzie and A. J. Roberts. Accurately model the Kuramoto- 
Sivashinsky dynamics with holistic discretisation. SI AM J. Applied 
Dynamical Systems, 5(3):365-402, 2006. 

[33] Tony MacKenzie. Create accurate numerical models of complex spatio- 
temporal dynamical systems with holistic discretisation. PhD thesis. 
University of Southern Queensland, 2005. 

[34] Tony MacKenzie and A. J. Roberts. Computer algebra derives the slow 
manifold of patch or element dynamics on lattices in two dimensions. 



Technical report. The University of Adelaide, 2011. http://arxiv.org/ 
liabs/1102.2037, 

[35] M. Marion and R. Temam. Nonlinear Galerkin methods. SIAM J. Nu- 
mer. Anal, 26(5):1139-1157, 1989. http : //locus . siam . org/SINUM/| 
[volume-26/art_0726063.htiiil , 



Tony Roberts, January 13, 2013 



References 



56 



[36] S. F. McCormick. Multilevel projection methods for partial differential 
equations. Society for Industrial and Applied Mathematics (SIAM), 
Philadelphia, 1992. 

[37] B. Mudavanhu and R. E. O'Malley. A new renormalization method for 
the asymptotic solution of weakly nonlinear vector systems. University 
of Washington, 2003. 

[38] National Physical Laboratory. Modern Computing Methods, volume 16 
of Notes on Applied Science. Her Majesty's Stationery Office, London, 
2nd edition, 1961. 

[39] G. A. Pavliotis and A. M. Stuart. Multiscale methods: averaging and 
homogenization, volume 53 of Texts in Applied Mathematics. Springer, 
2008. 

[40] W. H. Press, S. A. Teukolsky, W. T. Vetterhng, and B. P. Flannery. 
Numerical recipes in FORTRAN. The art of scientific computing. CUP, 
2nd edition, 1992. |http : // www. library . Cornell . edu/nr/cbookf pdf | 
fhtmlt 

[41] A. J. Roberts. Low-dimensional modelling of dynamics via computer 
algebra. Computer Phys. Comm., 100:215-230, 1997. .doi:10.1016/S0010] 
14655(96)00162-2 , 

[42] A. J. Roberts. Holistic discretisation ensures fidelity to Burgers' equation. 
Applied Numerical Modelling, 37:371-396, 2001. doi:10.1016/S0168-| 
||9274 (00) 00053-2, 

[43] A. J. Roberts. Derive boundary conditions for holistic discretisations of 
Burgers' equation. In K. Burrage and Roger B. Sidje, editors, Proc. of 
10th Computational Techniques and Applications Conference CTAC-2001, 
volume 44, pages C664-C686, March 2003. |http : // j ournal . a ustms | 



[org . au/ o j s/index . php/ANZIAMJ/ article/view/701 



[44] A. J. Roberts. A holistic finite difference approach models linear dynamics 
consistently. Mathematics of Computation, 72:247-262, 2003. 



Tony Roberts, January 13, 2013 



References 



57 



[45] A. J. Roberts. Holistic discretisation of dynamical partial differen- 



tial equations. Technical report, 


http : / / www . maths . adelaide . edu . au/ 


fanthony .roberts/holisticl .html 


April 2009. 



[46] A. J. Roberts. Holistic discretise three coupled dynamical partial differ- 
ential equations. Technical report, http : / /www . mat hs . adelaide . edu] 
jau/ anthony . roberts/holistic3 . html^ April 2009. 



[47] A. J. Roberts. Model dynamics across multiple length and time scales on 
a spatial multigrid. Multiscale Modeling and Simulation, 7(4): 1525-1548, 
2009. 

[48] A. J. Roberts and I. G. Kevrekidis. Higher order accuracy in the gap-tooth 
scheme for large-scale dynamics using microscopic simulators. In Rob 
May and A. J. Roberts, editors, Proc. of 12th Computational Techniques 
and Applications Conference CTAC-2004, volume 46 of ANZIAM J., 
pages C637-C657, July 2005. jhttp : / /anziamj . austms . org . au/V46/| 
[CTAC2004/Rob e [July 20, 2005]. 

[49] A. J. Roberts and I. G. Kevrekidis. General tooth boundary conditions for 
equation free modelling. SI AM J. Scientific Computing, 29(4):1495-1510, 
2007. 

[50] J. C. Robinson. The asymptotic completeness of inertial manifolds. 
Nonlinearity, 9:1325-1340, 1996. http : / / www . lop . org/E J/abstract/ 
10951-7715/9/5/013] 

[51] G. Samaey, I. G. Kevrekidis, and D. Roose. Damping factors for the gap- 
tooth scheme. In S. Attinger and P. Koumoutsakos, editors, Multiscale 
Modeling and Simulation, volume 39 of Lecture Notes in Computational 
Science and Engineering, pages 93-102. Springer- Verlag, 2004. 

[52] G. Samaey, I. G. Kevrekidis, and D. Roose. The gap-tooth scheme for 
homogenization problems. SIAM Multiscale Modeling and Simulation, 
4:278-306, 2005. ,doi:10.1137/030602046, 



Tony Roberts, January 13, 2013 



References 



58 



[53] R. Temam. Inertial manifolds. Mathematical Intelligencer, 12:68-74, 
1990. 

[54] H. A. van der Vorst. Iterative Krylov methods for large linear systems, 
volume 13 of Cambridge monographs on applied and computational math- 
ematics. Cambridge University Press, 2003. 

[55] A. Vanderbauwhede. Centre manifolds, normal forms, and elementary 
bifurcations. Dynamics Reported, 2:89-169, 1989. 

[56] S. Wolfram. The Mathematica Book, 3rd ed. Wolfram Media/ Cambridge 
University Press, 1996. 

[57] Keke Zhang and Gerald Schubert. Magnetohydrodynamics in rapidly 
rotating spherical systems. Annu. Rev. Fluid Mech., 32:409-443, 2000. 



Tony Roberts, January 13, 2013 



